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HIGHLIGHTS 


►  Time-dependent  mechanical  response  of  PFSA  membrane  is  investigated. 

►  Mechanical  response  of  PEM  fuel  cell  under  humidity  cycling  is  investigated  numerically. 

►  Time-dependent  model  predicts  higher  in-plane  stresses  compared  to  time-independent. 

►  Effect  of  humidity  cycle  variables  on  the  membrane  stresses  is  reported. 
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A  viscoelastic-plastic  constitutive  model  is  developed  to  characterize  the  time-dependent  mechanical 
response  of  perfluorosulphonic  acid  (PFSA)  membranes.  This  model  is  then  used  in  finite  element 
simulations  of  a  representative  fuel  cell  unit,  (consisting  of  electrodes,  gas  diffusion  layer  and  bipolar 
plates)  subjected  to  standardized  relative  humidity  (RH)  cycling  test  conditions.  The  effects  of  hold  times 
at  constant  RH,  the  feed  rate  of  humidified  air  and  sorption  rate  of  water  into  the  membrane  on  the  stress 
response  are  investigated.  While  the  longer  hold  times  at  high  and  low  humidity  lead  to  considerable 
redistribution  of  the  stresses,  the  lower  feed  and  sorption  rates  were  found  to  reduce  the  overall  stress 
levels  in  the  membrane.  The  redistribution  and  reduction  in  stress  magnitudes  along  with  inelastic 
deformation  during  hydration  eventually  lead  to  development  of  residual  tensile  stresses  after  dehy¬ 
dration.  Simulations  indicate  that  these  tensile  stresses  can  be  on  the  order  of  9—10  MPa  which  may  lead 
to  mechanical  degradation  of  the  membrane.  The  simulation  results  show  that  time-dependent  prop¬ 
erties  can  have  a  significant  effect  on  the  in-plane  stress  response  of  the  membrane. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Fuel  cells  represent  a  clean  alternative  to  current  technologies 
utilizing  hydrocarbon  fuel  resources.  However,  key  challenges  to 
wide-spread  commercialization  of  fuel  cell  technology  include  the 
high  cost  of  materials  and  the  relatively  low  durability  [1  ].  The  U.S. 
Department  of  Energy  (DoE)  has  set  a  durability  target  of  5000  h  of 
operation  for  PEM  fuel  cells  [2]  (corresponding  to  150,000  miles 
equivalent  in  a  vehicle).  Since  the  membrane  life  can  be  a  limiting 
factor  in  PEM  fuel  cell  durability,  this  work  focuses  on  the  durability 
of  the  PFSA  membranes  commonly  used  in  polymer  electrolyte 
membrane  (PEM)  fuel  cells. 

During  operation,  heat  and  water  are  generated  due  to  the 
electrochemical  reaction  associated  with  power  generation.  When 
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the  power  requirements  change  during  use,  the  fuel  cell  structure  is 
subjected  to  fluctuations  in  temperature  and  humidity.  As 
a  consequence,  the  structure  is  subjected  to  hygro-thermal  cyclic 
loading,  causing  stresses  in  the  membrane,  which  affect  the 
mechanical  durability. 

Chemical  and  mechanical  degradation  mechanisms  are  coupled 
in  actual  fuel  cell  operation.  Therefore,  a  full  understanding  of  the 
durability  of  PEM  fuel  cell  membranes  must  include  both  the 
chemical  and  mechanical  factors  that  contribute  to  material 
degradation  as  functions  of  time,  and  their  interaction  with  hygro- 
thermal  load  conditions  experienced  during  the  operating  and  non¬ 
operating  periods  [1].  However,  this  work  will  focus  on  the 
mechanical  aspects  of  degradation,  since  this  degradation  has  not 
been  studied  as  comprehensively  as  chemical  degradation. 

Mechanical  damage  is  typically  manifested  as  pinhole  formation 
or  cracking  and  can  be  accelerated  by  RH  cycling  as  well  as  low 
humidification  of  the  feed  stream  [3,4].  Once  these  flaws  are 
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formed,  gasses  may  pass  through  the  membrane  compromising  the 
fuel  cell  function  [3,4].  In  general,  the  mechanical  damage  is  caused 
by  the  cyclic  hygro-thermal  loading  in  the  cell  [5].  In  an  uncon¬ 
strained  state,  the  membrane  experiences  swelling  and  shrinkage 
during  water  sorption,  desorption  and  thermal  cycling.  However,  the 
membrane  electrode  assembly  (MEA)  is  constrained  from  defor¬ 
mation  by  the  stiffer  gas  diffusion  media  and  bipolar  plates  in  a  fuel 
cell.  These  constraints  cause  the  membrane  to  develop  compressive 
stress  at  higher  levels  of  water  uptake  and  high  temperature  during 
operation  [6-8].  If  these  compressive  stresses  are  high  enough  to 
cause  inelastic  deformation,  tensile  residual  stresses  may  develop 
upon  dehydration  and  cooling  of  the  membrane  [6,9]. 

There  are  several  cyclical,  hygro-thermal  accelerated  aging  tests, 
which  are  used  to  assess  the  mechanical  durability  of  membranes 
in  the  absence  of  electrochemical  activity  [2,10].  These  test  proto¬ 
cols  provide  a  basis  for  comparison  of  the  purely  mechanical 
behavior  of  different  membrane  material  compositions  and 
designs.  In  this  work,  we  will  focus  on  the  stress  response  of  PFSA- 
based  membranes  under  these  standardized  RH  cycling  conditions. 
In  order  to  determine  how  the  stresses  evolve  during  loading,  the 
mechanical  properties  of  the  PFSA  membrane  must  be  established. 

In  our  previous  work,  time-independent  linear-elastic-plastic 
models  were  used  to  characterize  the  long-term  membrane 
behavior  [6,9].  However,  our  recent  experiments  [11]  have  shown 
significant  time  dependence  that  may  have  a  significant  effect  on 
the  stresses  in  the  membrane  during  fuel  cell  loading.  Recent  work 
has  also  been  published  on  time-dependent  models  for  PFSA 
membranes  [12-14].  Most  of  the  models  assume  that  the  response 
is  viscoelastic,  therefore  inelastic  effects  have  been  not  been 
studied  thoroughly.  Patankar  et  al.  [15]  investigated  the  viscoelastic 
properties  of  Gore-Select®  57  PFSA  membrane  using  dynamic 
mechanical  analysis  (DMA).  They  used  time-temperature  super¬ 
position  to  create  a  hygro-thermal  master  curve.  Solasi  et  al.  [12] 
used  a  viscoplastic  constitutive  model  for  investigating  the 
behavior  of  a  dispersion  cast  PFSA  membrane.  They  measured  the 
uniaxial  stress  response  of  the  membrane  for  a  few  temperature 
and  hydration  levels  in  the  fuel  cell  operating  range.  The  model  did 
not  account  for  strain  rate  dependence  of  the  instantaneous  elastic 
modulus,  a  phenomenon  observed  in  their  experiments  and  in  our 
test  data  [11].  Silberstein  et  al.  [13]  used  a  micromechanically 
motivated  constitutive  model  to  capture  the  monotonic  and  cyclic 
stress  response  of  a  dispersion  cast  membrane.  The  measurements 
were  conducted  using  DMA  and  therefore  do  not  measure  the  long¬ 
term  stress-relaxation  behavior.  Yoon  et  al.  [14]  used  a  phenome¬ 
nological  constitutive  model  for  ionomer  membranes  based  on 
micromechanisms  of  polymer  deformation.  The  long-term  stress- 
relaxation  behavior  was  not  considered  while  validating  the  model. 
Yan  et  al.  [16]  characterized  the  material  response  of  the  PFSA 
membrane  using  a  non-linear  viscoelastic  model  based  on  Schap- 
ery’s  single  integral  method.  Lai  et  al.  [17]  used  a  linear  viscoelastic 
theory  with  time,  temperature  and  hydration  dependent  material 
properties  to  simulate  the  stresses  developed  in  the  membrane 
under  RH  cycling.  They  predicted  residual  tensile  stresses  in  the 
membrane  to  be  as  high  as  10.6  MPa,  which  is  near  or  higher  than 
the  proportional  limit  stress  at  various  conditions  [11]  suggesting 
the  need  for  a  viscoplastic  approach  to  modeling  of  the  membrane 
material.  Silberstein  et  al.  [18]  used  elastic- viscoplastic  properties 
for  the  PFSA  membrane  and  investigated  the  membrane  response 
via  experiments  under  partially  constrained  swelling  in  a  bimate¬ 
rial  strip  and  report  plastic  deformation  in  the  membrane.  Their  RH 
cycle  simulations  predicted  maximum  in-plane  stress  of  around 
12  MPa.  Li  et  al.  [19]  also  used  a  bimaterial  strip  to  measure  hygral 
stresses  in  the  membrane  and  attribute  the  tensile  stress  in  the 
membrane  upon  dehydration  to  the  viscoelastic  nature  of  the 
membrane  instead  of  plastic  deformation. 


From  the  above  discussed  studies,  it  appears  that  there  is  no 
consensus  regarding  plastic  deformation  in  the  membrane. 
However,  on  the  basis  of  our  experimental  data  [11  ]  and  numerical 
evaluation  of  strain  in  the  membrane  during  fuel  cell  RH  cycling 
(around  15%),  we  believe  that  there  is  inelastic  deformation  in  the 
membrane  and  it  can  play  a  critical  role  in  the  failure  mechanisms 
of  the  membrane.  Therefore,  a  viscoelastic-plastic  constitutive 
behavior  for  the  membrane  will  be  investigated  in  the  current 
work. 

In  a  companion  study  [11  ],  we  describe  the  set  of  uniaxial  tension 
and  stress-relaxation  tests  conducted  at  various  load  rates  and  hold 
strains  to  determine  the  viscoelastic-plastic  material  parameters 
over  a  range  of  environmental  conditions.  In  the  current  work,  we 
use  these  properties  to  develop  a  modified  two-layer  viscoelastic- 
plastic  constitutive  model  for  the  membrane.  Section  2  describes  the 
two-layer  viscoelastic-plastic  model  used  and  in  section  3,  the 
experimentally  determined  mechanical  properties  [11]  are  pre¬ 
sented  and  compared  to  the  constitutive  model  predictions.  Simu¬ 
lations  of  the  uniaxial  tensile  and  stress-relaxation  experiments  are 
discussed  in  Appendix  A  and  Appendix  B  respectively.  In  section  4 
we  incorporate  the  two-layer  model  into  a  finite  element  simula¬ 
tion  of  a  volume  element  in  a  typical  fuel  cell  to  explore  the  effects  of 
these  properties  on  the  stress  evolution.  Finally,  the  finite  element 
simulation  results  for  a  membrane  undergoing  two  standardized  RH 
cycle  tests  are  presented  and  discussed  in  section  5.  The  in-plane 
stresses  in  the  membrane,  which  have  been  shown  to  be  the  crit¬ 
ical  stresses  [7,8],  are  reported  for  a  selected  range  of  temperature 
and  humidity  cycle  time  parameters. 

2.  Constitutive  model  of  the  membrane 

In  this  section,  we  present  a  constitutive  model  capable  of 
capturing  the  hygro-thermo-mechanical  time-dependent  response 
of  the  membrane.  The  model  extends  the  linear-elastic  [9],  the 
linear-elastic,  perfectly-plastic  [6]  and  the  linear-elastic-plastic 
with  isotropic  hardening  [7]  models  we  have  previously  devel¬ 
oped  to  simulate  the  mechanical  response  of  the  membrane. 
The  constitutive  material  parameters  are  all  functions  of 
relative  humidity  (RH)  and  temperature  as  determined  in  the 
experiments  [11]. 

We  assume  that  the  total  strain  tensor,  ey,  is  given  as  the  sum 

ey  =  eV+el  +  ef  ;  (i  =  1,2,3),  (1) 

where  eh,  ejj  and  ejh  are  the  swelling,  thermal  and  mechanical  strain 
contributions,  respectively.  Each  of  these  strains  is  discussed  below. 


2.L  Swelling  strains 


Assuming  isotropic  swelling,  the  swelling  strains  in  the 
membrane,  caused  by  water  uptake,  are  calculated  using  the 
empirical  relationship  from  Kusoglu  et  al.  [20], 


(2) 


where  6  and  60  are  the  current  and  reference  temperature 
(expressed  in  degrees  K),  respectively.  </>P  is  the  polymer  volume 
fraction  given  by 


*  =  18A  T  EW/pp’  <3> 

where  A  is  the  number  of  water  molecules  attached  to  each  sul- 
phonic  acid  group,  EW  is  the  equivalent  weight  of  the  ionomer  and 
pP  is  the  density  of  the  dry  ionomer  material. 
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2.2.  Thermal  strains 


Assuming  isotropic  thermal  expansion,  the  thermal  strains  in  the 
membrane  caused  by  a  change  in  temperature  (0  -  0O)  are  given  by 


ejj  =  a(0  -  e0)5y,  (4) 

where  a  is  the  linear  coefficient  of  thermal  expansion  and  <5y  is  the 
Kronecker  delta  given  by 


5ij  = 


/ 1  4  =  i 

\o  ;l  *  J 


(5) 


2.3.  Mechanical  strains 


In  this  study,  the  mechanical  strains,  are  defined  in  the 
context  of  the  viscoelastic-plastic  model  [21,22]  shown  in  Fig.  1. 
This  model  consists  of  two  “layers”,  one  corresponding  to  the  time- 
independent  elastic-plastic  response  and  the  other  corresponding 
to  the  time-dependent  elastic-viscous  response.  The  two-layer 
model  was  chosen  due  to  its  relative  simplicity,  containing  6 
empirically  determined  material  parameters.  By  definition,  the 
total  mechanical  strain,  is  equal  in  both  layers 


M  _  , 

£ij  £ij 


.HP 


EV 


(6) 


where  the  superscripts  EP  and  EV  refer  to  elastic-plastic  and 
elastic-viscous  layers  of  the  model,  respectively.  In  the 
elastic-plastic  layer,  the  total  strain  is  the  sum  of  the  linear-elastic 
strain  in  the  spring  element  (efp)el  and  the  rate-independent  plastic 
strain  in  the  slider  elemental5  )pl 


EP 


(7) 


Extrapolating  the  spring  concept  in  the  elastic-plastic  layer  to 
represent  a  generalized  three-dimensional  isotropic  Hooke’s  law 
response,  the  elastic  strain  can  be  written  as  a  function  of  the  stress 
in  the  elastic-plastic  layer,  o|p 


1  +  vffEP  v  s  EP 

to' 


(8) 


Fig.  1.  One-dimensionai  idealization  of  the  two-layer  viscoelastic-plastic  model 
[21,22]. 


where  v  is  Poisson’s  ratio  and  JCP  is  the  stiffness  of  the  elastic 
element  in  this  layer. 

Assuming  isotropic  strain  hardening  (represented  by  the  spring, 
H,  parallel  to  the  slider  element),  the  yield  strength  is  a  function  of 
the  plastic  strain  as  well  as  being  a  function  of  the  relative 
humidity,  $1,  and  temperature,  6, 

=  ay((ej/P)pl’^’ (9) 

In  this  expression,  is  the  equivalent  plastic  strain  given  by 

(#)„,  -  /  no, 

where  d(^p)pl  are  the  increments  of  plastic  strain  tensor. 

Using  von  Mises’  yield  criterion  [23],  the  rate-independent 
plastic  flow  is  given  by 


(11) 


where  ay  is  the  yield  strength  of  the  material  and  Sy  are  the 
deviatoric  stress  tensor  components  defined  by 


%  —  r^kk^ij- 


(12) 


Under  von  Mises’  yield  criterion,  the  material  deforms  elasti¬ 
cally  for/((7y)  <  0  and  yielding  occurs  when 


>0.  (13) 

The  total  strain  in  the  elastic-viscous  layer  is  the  sum  of  linear- 
elastic  strain  in  the  spring  Kv  and  the  viscous  strain  in  the  dashpot 
element  given  by 


(14) 


where  subscript  el  refers  to  the  elastic  portion  and  vi  refers  to  the 
viscous  portion  of  the  strain.  Again,  assuming  isotropic  Hooke’s  law 
for  the  linear-elastic  constitutive  response  of  the  spring  in  this 
layer,  the  generalized,  three-dimensional  elastic  strain  as  a  function 
of  the  stress,  o|v,  is  given  by 


n-EV 


Ky 


V  r  1 

“  w6^ 


EV 
ijukk  • 


(15) 


For  the  dashpot  element,  the  Norton-Hoff  stress-strain  rate  law 
[24]  is  used,  the  generalized  three-dimensional  time-dependent 
form  of  which  is  written  as 


=  A(af 


where  A  and  n  are  material  parameters. 


(16) 


3.  Mechanical  properties  of  the  PFSA  membrane 

The  mechanical  response  of  PFSA  membrane  depends  on 
temperature,  water  content  and  loading  rate  [11].  In  this  section, 
the  experimental  data  from  the  uniaxial  tensile  and  stress- 
relaxation  tests  conducted  on  PFSA  Nafion®2111  membrane  are 
fitted  to  the  two-layer  model  introduced  in  the  previous  section 
using  the  equations  presented  in  Appendix  A  and  Appendix  B. 


1  Nafion  is  a  registered  trademark  of  E.I.  du  Pont  de  Nemours,  Inc.  Co. 
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The  material  parameters  thus  obtained  (Fig.  2)  were  used  in  FE 
model  simulations  of  the  test  specimens  to  compare  the  model 
predictions  to  test  data.  To  do  this,  a  plane  stress  ABAQUS  [22]  finite 
element  model  was  developed  with  gauge  length  of  50  mm  and 
width  of  20  mm,  and  meshed  using  CPS8R  elements  (plane  stress, 
8-node  biquadratic,  reduced  integration).  Selected  sample  results 
are  shown  in  Fig.  3.  The  figure  shows  that  the  constitutive  model 
captures  the  initial  elastic  and  non-linear  behavior  well.  We  assume 
that  the  onset  of  nonlinearity  corresponds  to  the  onset  of  inelastic 
behavior  of  the  membrane.  It  also  captures  the  variation  in  the 
initial  stiffness  and  the  overall  increase  in  stress  due  to  changing 
load  rates  fairly  well.  This  capability  is  not  present  in  the  original 
two-layer  model  [21,22],  but  was  accomplished  in  the  present  work 
by  adding  strain  rate  dependence  to  the  spring  element  in  the 
elastic-viscous  layer  of  the  viscoelastic-plastic  model.  The  strain 
rate  dependence  was  incorporated  into  the  model  by  introducing 
a  strain  rate  dependent  function  for  the  elastic-viscous  stiffness 
element,  I<v  and  using  values  for  the  stiffness  determined  from  the 
tensile  tests  conducted  at  various  load  rates.  In  the  finite  element 
model,  the  strain  rate  dependence  is  handled  via  an  ABAQUS  user 
subroutine  USDFLD  [22]  which  assigns  the  instantaneous  stiffness 
values  to  the  material  depending  on  the  current  strain  rate  in  the 
model. 

In  Fig.  4,  the  model  predictions  of  the  stress-relaxation  test  are 
compared  with  the  experimental  data  for  selected  conditions.  The 
model  captures  the  stress  peak  and  the  relaxation  behavior  for  the 
observed  time  frame  of  about  1—1.5  h.  The  trend  suggests  that 
relaxation  continues  beyond  the  test  time  although  at 
a  decreasing  rate.  Due  to  limitations  of  the  testing  equipment  and 
practical  considerations,  relaxation  beyond  1.5  h  was  not 
considered. 

Thus,  the  FE  simulations,  using  the  relatively-simple,  two-layer 
model  with  properties  shown  in  Fig.  2,  capture  the  experimentally 
measured  [11]  membrane  behavior  fairly  well.  In  the  next  section 
we  will  incorporate  these  material  properties  into  finite  element 
simulations  of  a  representative  volume  element  of  a  fuel  cell, 
subjected  to  humidity  cycling. 

4.  Numerical  implementation 

In  this  section,  the  effect  of  time-dependent  properties  on  the 
behavior  of  PFSA  membranes  in  a  fuel  cell  unit  operating  under  RFI 
cycling  is  investigated  via  numerical  simulations.  The  simulations 
are  based  on  finite  element  analyses  using  the  commercial  software 


Fig.  3.  [A]  True  stress  vs.  true  strain  uniaxial  tensile  response  of  the  PFSA  membrane  at 
various  environmental  conditions  obtained  experimentally  [2]  compared  with  the 
constitutive  model  predictions  for  strain  rate  3.3 -10-3  s-1  and  [B]  strain  rate 
8.3-10  2  s-1. 
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A 


Fig.  4.  [A]  True  stress  vs.  true  strain  stress-relaxation  response  of  the  PFSA  membrane 
at  various  environmental  conditions  obtained  experimentally  [2]  compared  with  the 
constitutive  model  predictions  for  hold  strain  0.05  and  [B]  hold  strain  0.1. 


ABAQUS  6.9  [22]  with  8-noded  biquadratic  coupled  temperature- 
displacement  generalized  plane  strain  elements  (CPEG8T)  and 
consists  of  23,805  elements  and  74,449  nodes  with  10  elements 
along  the  thickness  direction  of  the  membrane  (Fig.  5).  Numerical 
convergence  was  verified  using  progressively  finer  meshes  until  no 
appreciable  change  was  observed  in  the  results. 

4.1.  Geometry 

The  geometry  used  for  the  numerical  model  is  shown  in 
Fig.  5.  This  model  is  adapted  from  our  previous  work  [6,7,9] 
where  the  time-independent  response  was  studied.  General¬ 
ized  plane  strain  conditions  are  assumed,  corresponding  to 
a  uniform  strain  in  the  out-of-plane  (z)  direction.  A  typical  fuel 
cell  unit  contains  repetitive  grooves  and  lands  forming  the  gas 
channels  in  the  bipolar  plate.  We  constructed  a  model  for 
a  representative  repeating  element  consisting  of  a  half  of  a  land 
and  a  half  of  a  groove  with  the  dimensions  as  shown  in  Fig.  5. 


The  right  edge  and  bottom  have  symmetric  boundary  conditions 
imposed.  To  account  for  possible  displacement  in  the  in-plane 
(x)  direction,  a  continuous  boundary  condition  is  applied 
under  which  all  the  nodes  on  left  edge  of  the  model  can  move  in 
unison.  A  uniform  pressure  of  1  MPa  is  applied  from  the  top 
edge  of  bipolar  plate  to  simulate  a  typical  spring-loaded 
clamping  pressure  for  a  fuel  cell. 

4.2.  Properties 

The  following  material  properties  and  behaviors  are  used  for  the 
various  components  in  the  finite  element  model: 

1.  Viscoelastic-plastic  properties  derived  from  the  experiments 
are  used  for  the  PFSA  membrane  (section  3).  The  material 
properties  and  the  sorption  coefficient  for  the  membrane  are 
incorporated  into  the  model  as  functions  of  water  content  and 
temperature.  Poisson’s  ratio  of  0.4  is  assumed  for  the 
membrane  [25]. 

2.  The  membrane  is  the  only  component  that  is  assumed  to  swell 
due  to  water  absorption. 

3.  Thermal  expansion  and  swelling  of  the  membrane  are  assumed 
isotropic  for  simplicity. 

4.  The  bipolar  plates  are  assumed  to  be  linear-elastic  and  isotropic 
with  elastic  modulus  and  Poisson’s  ratio  of  10  GPa  and  0.25 
respectively.  Thermal  expansion  is  neglected. 

5.  The  gas  diffusion  layer  (GDL)  is  assumed  to  be  linear-elastic  and 
transversely  isotropic  with  in-plane  and  the  out-of-plane 
elastic  moduli  of  1500  MPa  and  9  MPa  respectively  [26].  The 
coefficient  of  thermal  expansion  is  2.6  x  1CT9  K-1. 

6.  The  electrodes  are  assumed  to  have  time-independent 
isotropic,  elastic-plastic  response  with  isotropic  hardening. 


4.3.  Load  cases 

Standard  test  protocols  have  been  established  by  the  US 
Department  of  Energy  (DoE)  to  assess  the  performance  and 
durability  of  fuel  cell  components  under  cyclic  loading  [2,27].  The 
test  protocol  that  targets  mechanical  durability  is  used  for  inves¬ 
tigation  of  membrane  performance  in  this  work.  This  cycle  is 
referred  to  as  the  “DoE  cycle”  in  Fig.  6  and  involves  a  constant 
temperature  of  80  °C  while  the  humidity  is  cycled  equally  on  both 
the  anode  and  cathode  sides,  between  dry  and  supersaturated 
conditions  (corresponding  to  the  dewpoint  at  90  °C).  The 
maximum  water  volume  fraction  used  in  the  analysis  is  O  =  0.32. 
Using  Eq.  (3),  this  corresponds  to  a  water  content  in  the 
membrane  (7)  of  around  14.  The  cell  is  held  at  the  dry  and  wet 
state  for  2min  each,  and  the  transition  is  assumed  to  take  place  in 
a  very  short  time  (Is  in  the  simulation). 

An  alternative  RH  cycling  test,  developed  by  W.  L.  Gore, 
referred  to  as  the  “Gore  cycle”  in  Fig.  6,  is  also  investigated  in  this 
work  [28].  This  cycle  involves  similar  conditions  as  the  DoE  cycle 
except  that  the  cell  is  held  in  the  wet  state  for  10s  and  in  the  dry 
state  for  50s. 

In  both  tests,  humid  and  dry  air  are  alternately  forced  into  the 
flow  channels,  and  therefore  the  actual  water  profile  at  the 
membrane  is  mediated  by  the  diffusion  through  GDL,  electrode  and 
the  membrane  itself.  The  water  volume  fraction  in  the  membrane 
and  the  other  fuel  cell  components  is  modeled  as  an  additional 
degree  of  freedom  in  addition  to  the  nodal  displacements  in  the 
finite  element  analysis.  The  water  volume  fraction  for  the 
membrane  at  the  interface  with  the  electrodes  has  a  time-lag  as 
shown  in  Fig.  7.  Snapshots  of  the  water  content  contours  during 
hydration  and  dehydration  are  shown  in  Fig.  8. 
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Fig.  5.  Schematic  of  the  2D  generalization  of  fuel  cell  unit  (not  to  scale). 


5.  Results  and  discussion 

In  our  previous  work  [6,20,29]  we  have  shown  that  the  in-plane 
stress  is  the  dominant  stress  in  the  membrane  during  fuel  cell 
operation.  Consequently,  we  will  focus  on  that  component  only.  The 
in-plane  stress  varies  both  as  a  function  of  time  and  location.  The 
snapshots  displaying  the  spatial  variation  of  the  in-plane  stresses  in 
the  membrane  during  the  fuel  cell  simulations  under  standard  RH 
cycle  loading  condition  (section  4.3)  are  shown  in  Fig.  9.  The 
evolution  of  the  in-plane  stress  in  the  membrane  for  an  element  on 
the  groove  side  of  the  fuel  cell  is  presented  in  Fig.  10.  The  particular 
‘observation  point’  as  chosen  in  Fig.  10  is  a  representative  element 
selected  to  show  the  general  stress  evolution  behavior.  The  in-plane 
stresses  caused  by  the  initial  clamping  pressure  of  1  MPa,  applied  at 
the  top  of  bipolar  plate  are  small  in  magnitude  (Fig.  9A).  A  parametric 
study  was  conducted  to  investigate  the  effect  of  clamping  pressure 


(not  shown  for  brevity).  The  increased  pressure  up  to  3  MPa  results 
in  a  minor  increase  (<5%)  in  the  compressive  stresses  after  clamping 
and  hydration  and  therefore,  a  similar  decrease  was  found  in  the 
residual  tensile  stress  after  dehydration. 

When  the  system  is  subjected  to  hydration,  the  membrane  takes 
on  water  and  attempts  to  swell  (from  A-B  in  Fig.  10).  Due  to  the 
mechanical  constraints,  the  membrane  is  mostly  prevented  from 
deforming  in  the  in-plane  direction,  resulting  in  relatively  large 
compressive  stresses,  on  the  order  of  9-10  MPa  (Fig.  9B).  During 
hold  at  high  humidity  (10s  duration  from  point  B  to  C  in  Fig.  10), 
these  compressive,  in-plane  stresses  begin  to  relax.  The  redis¬ 
tributed  stress  profile  is  shown  in  Fig.  9C  and  it  is  clear  that  the 
overall  level  of  compressive  stress  is  reduced,  along  with  a  signifi¬ 
cant  redistribution  of  the  stresses.  After  hold  at  high  humidity,  the 
system  is  flushed  with  dry  nitrogen  (C-D  in  Fig.  10)  resulting  in 
a  decrease  in  water  content  and  a  decrease  in  compressive  stresses. 
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Time  [s] 

Fig.  6.  DoE  and  GORE  RH  Cycling  Test  profile. 


This  decrease  eventually  leads  to  residual  tensile  in-plane  stresses, 
on  the  order  of  7— 8  MPa,  as  shown  in  Fig.  9D.  During  the  holding  at 
dry  conditions  (50s),  the  membrane  in-plane  tensile  stresses  relax 
considerably  (D— E  in  Fig.  10)  and  are  finally  redistributed  as  shown 
in  Fig.  9E.  The  figures  show  that  residual  stresses  under  the  groove 
are  higher  than  those  under  the  land.  This  confirms  the  trend  from 
our  previous  work  with  time-independent  properties  [6,9,20,29]. 
However,  the  time-dependent  properties  implemented  here,  result 
in  larger  peak  for  the  in-plane  stresses  after  hydration  and  dehy¬ 
dration  compared  to  the  time-independent  predictions.  During  the 
hold  phases,  these  peak  stresses  relax  to  levels  below  those  seen  in 
time-independent  simulations. 

The  hold  time  at  constant  hydration  and  dehydration  is  longer  in 
the  DoE  RH  cycle  than  the  Gore  RH  cycle,  allowing  for  more 
relaxation  (Fig.  10).  This  results  in  larger  residual  tensile  stresses 
after  dehydration  at  the  end  of  the  cycle.  However,  the  Gore  RH 
cycling  test  involves  a  higher  frequency  of  loading,  which  produces 
more  reversals  of  the  stresses  for  a  given  period  of  testing.  There¬ 
fore,  without  more  information  on  the  fatigue  resistance  of  the 
membranes,  it  is  not  possible  to  determine  from  these  simulations 
which  test  cycle  will  cause  earlier  failure  in  the  membrane. 

The  typical  strain  evolution  in  the  membrane  during  Gore  RH 
cycling,  for  a  representative  element  on  the  groove  side  is  shown  in 

Fig.  11.  The  plastic  strain  magnitude,  ep  =  yj 2/3^  •  eP1  is  about 
0.15,  and  therefore,  can  be  instrumental  in  eventual  failure  of  the 
membrane. 

The  transition  time  (that  is,  the  time  assumed  for  the  input  RH 
load  to  change  linearly  from  dry  to  saturated  conditions  in  the 
standardized  RH  cycles)  is  not  clearly  defined  by  the  testing 
protocol  and  may  vary  for  different  regions  of  the  fuel  cell  stack. 
Therefore,  we  will  study  the  effect  of  this  transition  by  varying  this 
time  and  observing  the  resulting  in-plane  stress  evolution.  The 
peak  in-plane  compressive  and  tensile  stresses  (corresponding  to 
points  B  and  D  in  Fig.  10,  respectively)  for  the  cycle  simulations 
conducted  using  selected  transition  times  are  shown  in  Fig.  12.  The 
results  indicate  that  longer  transition  times  result  in  reductions  of 
both  the  peak  compressive  stresses  after  hydration  and  peak  tensile 
stresses  after  dehydration. 

In  the  standard  RH  cycle  protocols,  the  humidified  air  is  fed 
simultaneously  into  both  the  anode  and  the  cathode  side  of  the  fuel 
cell.  However,  during  real  fuel  cell  operation,  the  anode  side  is  often 
less  hydrated  due  to  the  production  of  water  at  the  cathode  and  the 
electro-osmotic  drag,  creating  a  water  content  gradient  [30].  To 
investigate  the  effect  of  loading  with  gradient  water  content,  we 


Fig.  7.  Water  volume  fraction  variation  as  seen  at  membrane-electrode  interface 
compared  to  input  at  flow  channels  indicating  swelling  time-lag. 
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Fig.  8.  FE  results  showing  water  content  in  the  MEA  at  the  end  of  [A]  hydration  and  [B] 
dehydration  step  of  the  RH  cycle. 
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Fig.  9.  Contours  of  the  in-plane  stress  in  the  membrane  during  Gore  RH  cycling  (first 
cycle);  [A-E]  refer  to  time  instants  as  in  Fig.  10. 


simulated  the  case  where  the  cathode  side  of  the  membrane  is 
continually  subjected  to  saturated  air,  while  the  air  at  the  anode 
side  is  cycled  according  to  the  Gore  RH  cycle  at  a  temperature  of 
80  °C.  The  contours  displaying  the  spatial  variation  of  the  in-plane 
stress  through  the  thickness  are  given  in  Fig.  13.  Asymmetric 
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Fig.  11.  Typical  strain  evolution  during  Gore  RH  cycling  for  an  element  in  the 
membrane  on  the  groove  side.  In-plane  swelling,  elastic,  viscous  and  plastic  strain  is 
shown  along  with  total  plastic  strain  magnitude. 


distribution  of  in-plane  stress  due  to  non-uniform  water  content  is 
observed  (Fig.  13B,  C).  After  dehydration,  a  mixed  state  of  stress  is 
observed  with  the  portion  of  the  membrane  under  the  land  in 
compression  and  the  portion  under  the  groove  experiencing 
tension  (Fig.  13C). 

To  further  explore  this  phenomenon,  the  water  content  varia¬ 
tion,  in-plane  strain  and  resulting  in-plane  stress  evolution  for  two 
elements  under  the  groove,  one  on  the  anode  and  one  on  the 
cathode  side  of  membrane  are  shown  in  Fig.  14.  It  is  evident  that 
due  to  the  continual  hydration  at  the  cathode,  the  membrane 
experiences  a  smaller  differential  of  water  content  and  less 
swelling-shrinkage  strain,  leading  to  lower  overall  stress. 

As  the  humidified  and  dry  air  are  alternatively  fed  into  the 
channels,  the  sorption  process  limits  the  rate  of  change  in  the  water 
content  in  the  GDL,  catalyst  and  membrane.  In  the  present  work, 
the  sorption  phenomenon  is  modeled  using  a  diffusion  coefficient 
for  the  membrane  as  an  empirically  derived  temperature  and  water 


Time  [s]  Time  [s] 


Fig.  10.  In-plane  stress  evolution  comparing  elastic-plastic  to  viscoelastic-plastic  properties  during  [A]  Gore  and  [B]  DoE  RH  cycling  test  in  a  groove  side  element  of  the  membrane. 
The  markers  [A-E]  correspond  to  the  stress  contours  shown  in  Fig.  9. 
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Fig.  12.  Effect  of  hydration/dehydration  time  on  the  in-plane  compressive  and  tensile 
stresses  in  the  membrane  for  Gore  and  DoE  RH  cycling. 


content  dependent  function  taken  from  the  work  of  Weber  and 
Newman  [31  ].  This  coefficient  affects  the  dynamic  water  uptake  for 
the  membrane  that  in  turn  affects  the  stress  levels.  In  order  to 
investigate  this  effect,  we  ran  the  simulation  with  different  frac¬ 
tional  values  of  this  coefficient  [31  ].  Since  the  original  coefficient 
values  from  the  literature  [31]  give  relatively  fast  water  uptake, 
smaller  (slower)  coefficient  values  were  chosen  for  the  parametric 
analysis.  The  resulting  water  profile  variation  is  shown  in  Fig.  15A 
and  the  in-plane  stress  evolution  in  Fig.  15B.  The  results  indicate 
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Fig.  13.  In-piane  stress  contours  in  the  membrane  under  non-uniform  (gradient) 
loading.  [A-E]  refer  to  time  instants  in  Fig.  14. 
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Fig.  14.  Water  volume  fraction,  in-plane  strain  and  in-plane  stress  variation  for  an 
element  under  the  groove  on  the  anode  side  and  one  on  the  cathode  side  of  the  PFSA 
membrane  for  the  first  cycle  under  non-uniform  (gradient)  loading  the  markers  [A-E] 
correspond  to  the  stress  contours  shown  in  Fig.  13. 


Time  [s] 

Fig.  15.  For  three  values  of  sorption  coefficient  D(4>,T)  obtained  from  Weber  et  al.  [31  ], 
[A]  Water  volume  fraction  variation  in  the  membrane  during  RH  cycling;  [B]  In-plane 
stress  evolution. 
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that  with  a  decreasing  value  of  the  coefficient,  D,  the  membrane 
never  gets  fully  saturated  (corresponding  to  water  volume  fraction, 
0  =  0.32)  during  the  hydration  phase  of  the  cycle  and  therefore,  the 
stress  levels  are  significantly  reduced.  In  addition,  for  the  lower 
sorption  rate  D/50,  the  membrane  does  not  reach  the  dry  state  ($ 
=  0.06)  during  the  dehydration  phase,  further  reducing  the  stress 
levels. 

6.  Concluding  remarks 

In  this  work,  a  two-layer  viscoelastic-plastic  constitutive  model 
[21,22]  was  adapted  to  characterize  the  time-dependent  and  rate- 
dependent  mechanical  response  of  PFSA  membrane  under  various 
environmental  conditions  and  load  rates  spanning  a  range  of  fuel 
cell  operation  conditions.  We  modified  the  standard  constitutive 
model  [21,22]  by  adding  strain-rate  dependence  and  determined 
the  specific  time-temperature-humidity  dependent  parameters 
from  our  tensile  and  stress-relaxation  experimental  data  [11]  for 
Nafion®211  membrane.  The  simulations  in  section  2  show  that  the 
resulting  model  captures  the  material  response  well  over  the  range 
of  temperature,  humidity  and  load  rates  studied. 

The  constitutive  model  was  incorporated  into  a  finite  element 
analysis  to  study  the  stresses  developed  in  a  PFSA  membrane, 
including  time-dependent  effects,  during  accelerated  mechanical 
testing.  The  analysis  was  conducted  as  a  coupled  sorption- 
displacement  problem  [20]  enabling  us  to  model  the  water 
content  profile  and  mechanical  stresses  simultaneously. 

The  results  of  the  simulations  show  that  (i)  the  time-dependent 
models  predict  larger  peak  compressive  and  residual  tensile 
stresses  in  the  membrane  as  compared  to  earlier  predictions  using 
time-independent  properties;  (ii)  slower  feed  rates  for  the 
humidified  air  reduce  the  overall  in-plane  stresses;  (iii)  slower 
sorption  rates  decrease  the  swelling  strain  considerably  and  hence 
lower  the  stresses;  (iv)  residual  in-plane  tensile  stresses  developed 
in  the  DoE  RH  cycles  are  higher  than  for  the  Gore  RFI  cycle,  owing  to 
the  longer  relaxation  at  constant  strains  in  the  former;  (v)  Uniform 
moisture  feed  on  both  sides  of  the  membrane  produces  higher 
stresses  as  compared  to  non-uniform  loading. 

In  conclusion,  these  results  suggest  that  the  mechanical  degra¬ 
dation  and  hence  the  membrane  lifetime  may  be  influenced  by  the 
inherent  time  rates  of  humidity  variation  during  real  fuel  cell 
operation.  The  sites  of  the  largest  stress  magnitude  and  cyclic 
reversals  may  act  as  potential  damage  initiation  points  and  trigger 
ultimate  failure  of  the  membrane.  This  study  adds  to  the  under¬ 
standing  of  the  mechanical  aspects  of  durability  of  PEM  fuel  cells. 
The  numerical  model  along  with  carefully  determined  material 
properties  helps  ascertain  severe  loading  cases  for  the  hygro- 
thermal  cycling  in  a  fuel  cell. 
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Appendix  A.  Tensile  simulation 

In  this  section,  derivation  of  the  stresses  in  individual  layers  of 
the  two-layer  viscoelastic-plastic  model  [21,22]  for  tensile  loading 
is  explained.  In  order  to  visualize  the  individual  effects  of  the 
parameters  in  the  two-layer  model  (see  Fig.  1),  we  investigate  the 
model  response  during  a  generic  tensile  test  simulation,  using 
generic  properties. 

In  a  displacement-controlled  tensile  test,  the  strain  rate,  e,  is 
nearly  constant.  The  strain  evolution  for  each  of  the  individual 
parameters  of  the  constitutive  model  is  shown  in  Fig.  Al-A.  In  the 


elastic-plastic  layer,  the  strain,  eEf,  in  the  elastic  element,  I<P, 
changes  linearly  until  the  elastic  limit  is  reached.  Thereafter,  the 
elastic  strain  increases  monotonically  at  a  reduced  rate.  The  plastic 
strain,  £Ej\  appears  after  the  proportional  limit  is  reached  and 
increases  linearly  beyond  this  limit.  Similarly,  in  the  elastic-viscous 
layer,  the  strain,  e EjV,  in  the  elastic  element,  JCV,  initially  increases 
linearly  until  the  time-dependent  dashpot  is  activated.  Once  the 
dashpot  is  activated,  there  is  no  further  increase  in  the  strain  for  the 
elastic  element,  while  the  viscous  strain,  eEV,  increases  linearly.  For 
one-dimensional  loading,  the  total  stress  is  given  by  the  sum  of 
stresses  in  the  two  layers  (Fig.  Al-B), 


The  stress  (jEP  in  the  elastic-plastic  layer  is  given  by 


„EP  _  /  Kpe  ;  aEP  <  <7y 

l(7y+^(e  —  £c)  ;<7EP>(7y 


A-2 


where  ec  =  oy/Kp,  /<p  is  the  elastic  stiffness  and  H  is  the  post-yield 
slope  of  the  true  stress  vs.  true  strain  curve  for  the  elastic-plastic 
layer.  H  can  generally  be  any  function  of  strain  in  addition  to 
having  hygral  and  thermal  dependence.  Using  the  information  from 
our  experiments,  the  form  of  H  used  in  this  work  is  such  that 


H  =  H(e,  A,  6)  =  H(A,  0)  ;  e  <  0.3.  A-3 

The  stress  in  the  elastic-viscous  layer  results  in  a  first-order 
inhomogeneous,  non-linear  differential  equation.  Again,  assuming 
one-dimensional  loading,  the  individual  elastic  and  viscous  stresses 
in  the  elastic-viscous  layer  are  equal  to  the  total  stress 


tev 


EV 


(J  —  <7„i  —  (T, 


-EV 
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The  strain  in  the  elastic- viscous  layer,  given  by  Eq.  (14)  when 
written  in  one-dimensional  form  (ij  =  1)  gives  the  following 
expression 


EV 


_  pEV 
£el  ' 


EV 
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where  the  elastic  strain  is  related  to  the  stiffness,  I<v  and  is  given  by 


aEV 
CEV  _  ael 


I<V  ' 


A- 6 


Differentiating  Eq.  A-5  with  respect  to  time,  t,  and  substituting 
the  viscous  strain  eEV  from  Eq.  (16)  and  elastic  strain  from  Eq. 
(A-6),  we  get  a  relationship  between  the  strain  rate,  stress  and 
stress  rate, 


•  EV 

*EV  =  ^L  +  /1 


Kv 


A- 7 


where  a  “dot”  above  the  quantity  represents  its  derivative  with 
respect  to  time.  Rearranging  Eq.  (A-7),  and  using  Eq.  (A-4)  we  can 
write  the  resulting  constitutive  equation  for  the  elastic-viscous 
layer  as 


<7EV+/\i<:v(ffEV)n-/<v£EV  =  o, 


A-8 


where  eEV  =  eM  =  e  is  the  total  strain  rate.  The  solution  to  this 
equation  is  found  numerically  using  illustrative  properties  and  the 
resulting  stress  (Fig.  Al-B)  shows  an  increase  with  an  increase  in  I<y 
to  a  constant  maximum  value  of 


l/n 


A-9 
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Eq.  (A-9),  when  combined  with  Eq.  (A-2)  gives  the  total  stress 
for  the  model  under  constant  strain  rate  loading  (assuming 
inelastic  deformation)  as 

/A  1/11 

a  =  <ry  +  H(e  -  ec)  +  ( j  J  .  A-10 


Appendix  B.  Relaxation  simulation 

In  a  typical  stress-relaxation  test,  the  specimen  is  quickly  loaded 
to  a  desired  strain.  The  strain  is  then  held  constant  and  the 
reduction  in  stress  is  measured  as  a  function  of  time.  During 
a  stress-relaxation  simulation  using  generic  properties,  the  strain  in 
the  individual  components  of  the  constitutive  model  varies  as 
shown  in  Fig.  Bl-A.  The  stress  in  the  model  at  the  end  of  the  initial 
loading  (tension)  can  be  determined  by  Eq.  (A-10),  which  is  also  the 
initial  stress  for  the  relaxation  phase.  During  the  relaxation  period, 
the  stress  in  the  elastic-plastic  layer,  (jEP,  is  unchanged,  while  the 
stress  in  the  elastic-viscous  layer,  <jev,  decreases  due  to  dissipation 
in  the  viscous  element,  as  illustrated  in  Fig.  Bl-B.  The  stress  in  the 
elastic-plastic  layer  is  still  given  by  Eq.  (A-2),  but  during  relaxation, 
the  stress  in  the  elastic-viscous  layer  evolves  as  the  strain  compo¬ 
nents  are  redistributed.  Therefore,  to  determine  total  stress  for  the 
model,  the  stress  in  the  elastic-viscous  layer  alone  needs  to  be 
evaluated. 

During  the  stress-relaxation  period,  the  strain,  eE^,  in  the  elastic 
element,  Kv,  of  the  elastic-viscous  layer  decreases  as  the  dashpot 
lengthens,  such  that  the  total  strain  in  the  layer  remains  constant 
(Fig.  Bl-A), 


_EV  ,  E V  _  p 
£e\  +£vi  -  £- 


B-l 


Thus,  upon  differentiating  and  rearranging  Eq.  (B-l),  the 
following  relationship  between  the  strain  rates  is  obtained, 


•  EV 


_  aEV 

~£e\  ■ 


B-2 


recovered  from  the  Eq.  (B-6)  by  substituting  t  =  0.  During 
stress-relaxation,  which  occurs  for  t  >  0,  the  stress  in  the  elastic- 
viscous  layer,  aEV,  given  by  Eq.  (B-5)  decreases  as  shown  in 
Fig.  Bl-B,  and  as  a  result  the  total  stress,  a  given  by  Eq.  (B-6) 
relaxes  to  an  equilibrium  value  equal  to  the  stress  in 
elastic-plastic  layer  given  by  Eq.  (A-2). 

The  stress-relaxation  data  obtained  at  each  environmental 
condition  is  fit  to  the  Eq.  (B-6)  and  the  elastic-viscous  parameters  A, 
n  and  I<v  are  determined.  The  time-independent  portion  of  the  Eq. 
(B-6)  which  is  also  given  by  Eq.  (A-2),  is  then  used  to  construct 
a  ‘quasi-static’  curve  from  which  the  elastic-plastic  parameters  oy, 
H  and  I<p  are  determined. 


The  viscous  strain  rate  (e^)  is  governed  by  Eq.  (16),  in  which  it 
can  be  replaced  by  the  corresponding  negative  elastic  strain  rate 
using  Eq.  (B-2),  giving 

-e™  =  A^)".  B-3 

Eq.  (B-3)  can  be  simplified  using  Eqs.  (A-4)  and  (A-6)  to 
ffEV+A/<v((TEV)n=  0.  B-4 


Eq.  (B-4)  is  a  first-order  homogeneous  non-linear  equation  in 
(jEV  with  the  solution 


<7EV  = 


n  -  1  \AKyt 


B-5 


Therefore  the  total  stress  during  stress-relaxation  is  given  by  the 
sum  of  Eqs.  (A-2)  and  (B-5)  as 


(7  —  (7y  H(^£  —  £“c) 


(1-n) 


n  -  1  )AKvt 


B-6 


The  maximum  value  of  the  total  stress  a  at  the  end  of 
a  tensile  pull  with  strain  e  as  given  by  Eq.  (A-10)  can  be 


Fig.  Al.  [A]  Typical  strain  evolution  and  [B]  typical  stress  evolution  during  uniaxial 
displacement-controlled  tensile  test  simulation. 
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Fig.  Bl.  [A]  Typical  strain  evolution  and  [B]  typical  stress  evolution  during  uniaxial 
stress-relaxation  test  simulation. 
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